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I. INTRODUCTION 

The progressive miniaturization of electronic devices 
and their future application in nanoscale circuitry entails 
the necessity of developing a quantum theory of trans- 
port. Such a theory should account for the full atomistic 
structure of the contacts and the molecular device, and 
for the dynamical effects of exchange and correlation on 
the electron motion. 

Most theoretical works on quantum transport focus 
on steady-state properties and are based on a "con- 
tacting approach" first introduced by Caroli et al. In 
this approach the system is separated in isolated parts 
in the remote past (left /right leads and central device) 
and the contacts between subsystems are treated as a 
time- dependent perturbation (adiabatic switching). Such 
a partition is crucial in the formal derivation of the 
Meir-Wingreen formula, 5 which allows for studying time- 
dependent phenomena and correlation effects in the de- 
vice region but is not suitable to include long-range in- 
teractions and memory effects in a first-principle manner 
(see discussion below) . 

The contacting approach has a severe limitation. In a 
typical experiment the whole system is contacted and in 
equilibrium before a driving external field is switched on. 
Thus, there is an implicit assumption of equivalence be- 
tween a) the initially partitioned and biased system once 
the contacts are introduced and b) the whole system in 
equilibrium once the bias is turned on. This assumption 
might be reasonable in the long-time limit but is clearly 
inappropriate at any finite time. As a consequence, the 
contacting approach is not suitable for describing tran- 
sient regimes. Also, the memory of initial conditions is 
not properly accounted for. In a non-interacting system 
memory effects are washed out provided the local density 
of states is smooth in the device region. 6-8 However, this 
is generally not true in the interacting case. 7 ' 9,10 

Recently Dhar and Sen 11 have shown that memory ef- 
fects might be observed also in non-interacting systems 
with bound states. This finding renders the contacting 
approach ambiguous since the equilibrium distribution of 



the initially isolated device is arbitrary. It has also been 
shown that the one-particle density matrix of the unbi- 
ased system does not reduce to the standard equilibrium 
result but rather oscillates with frequencies lu = Sb — Ew , 
where Eb, Sb> are bound-state energies. In order to cir- 
cumvent this problem Dhar and Sen proposed to intro- 
duce two "extra" reservoirs weakly coupled to the device. 
For the unbiased system the equilibrium result is recov- 
ered in the limit where the "extra" couplings tend to 
zero. However, out of equilibrium the results depend on 
the limiting procedure. 

In this work we use the "partition-free approach" 
by Cini 12 which is free from all above limitations. 
In contrast to the contacting approach the system is 
not partitioned in the remote past and is in equi- 
librium at a unique temperature and chemical poten- 
tial (thermodynamic consistency). The initial equilib- 
rium distribution of the device is unambiguous. The 
system is driven out of equilibrium by exposing the 
electrons to a time- dependent electric field. Thus, 
the external perturbation is a local potential and the 
partition-free approach can be combined with Time- 
Dependent Density Functional Theory 13-16 (TDDFT) 
to calculate total currents and densities in interacting 
systems. 7,8 The use of TDDFT in quantum transport is 
gaining ground 7,8, 10:17-29 and several properties of the 
time-dependent exchange-correlation potential and ker- 
nel have recently been discussed. 30-34 

In a previous work 6 we have shown how a steady cur- 
rent develops under the influence of a constant bias. For 
non-interacting electrons we also proved that the steady 
current is independent of the history of the applied bias 
and agrees with the steady current calculated in the con- 
tacting approach. The theory has been developed as- 
suming a smooth density of states in the device region. 
Here, we generalize the theory and include bound states 
in the description of time-dependent quantum-transport 
phenomena. Our main findings are: For non-interacting 
electrons 1) in the presence of bound states the total cur- 
rent and the one-particle density matrix oscillate with 
frequencies uj = Eb — £&', and the amplitude of the os- 
cillations are unambiguous calculable quantities, 2) the 
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amplitude of the oscillations depends on the history of 
the applied bias and on the original equilibrium, mean- 
ing that the long-time limit does not wash out the effects 
of different initial conditions, 3) for the unbiased sys- 
tem the oscillations vanish and all the time-dependent 
quantities reduce to their equilibrium value. Thus, in 
the partition-free approach there is no need of "extra" 
reservoirs to recover standard equilibrium properties. For 
Kohn-Sham electrons (TDDFT) 4) the steady-state as- 
sumption is not consistent with the presence of bound 
states, 5) it is possible to have self-consistent oscillatory 
solutions even for time- independent external fields, 6) 
density oscillations might induce oscillations in the ef- 
fective potential of TDDFT and hence give rise to new 
conductive channels. This effect can not be captured in 
static DFT calculations. 

The plan of the paper is as follows. In Section II we 
give a brief introduction to the extended Keldysh for- 
malism which properly accounts for the effects of initial 
correlations. A general formula for the total current is 
derived in Section III. Such a formula explicitly contains 
the distribution function of the system in equilibrium and 
allows for working in both the contacting approach and 
the partition-free approach. In Section IV we generalize 
the theory of Ref. 6 to include bound states in quantum 
transport. The implications of our findings in a formu- 
lation based on TDDFT are discussed in Section V. In 
Section VI we summarize the main results and draw our 
conclusions. 



II. THE KELDYSH-GREEN'S FUNCTION 

Let us consider a system of non-interacting electrons 
(or mean-field electrons or Kohn-Sham electrons) de- 
scribed by 
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JJ(t)=£[Jf(i)] 



(i) 



where the matrix 



H(t) = 



v CL (t) 





v LC (t) 
£c(t) 
v RC {t) 





v CR (t) 
£ R (t) 



(2) 



models the electrode-device-electrode system of Fig. 1. 
(Here and in the following we use boldface to indicate 
matrices in one-electron labels.) 

Without loss of generality we assume that the system 
is in equilibrium for negative times t < 0. Letting fi be 
the chemical potential and (3 be the inverse temperature, 
all equilibrium quantities can be expressed in terms of 
the density matrix 



P : 



(3) 



with N = J2m c m c m the operator of the total number of 
particles and H° = H(t < 0). The Matsubara-Green's 
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FIG. 1: The electrode-device-electrode system described by 
H(t). It consists of left (L) and right (R) metallic regions 
coupled to a central scattering region C. The region C con- 
tains the actual molecular device and few atomic layers of the 
L/R electrodes. 



function technique is a well established theory for cal- 
culating /5-averaged quantities but is limited to equilib- 
rium problems. A very powerful tool for dealing with 
non-equilibrium (as well as equilibrium) problems is pro- 
vided by the non-equilibrium Green's function (NEGF) 
theory. 35 ' 36 
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FIG. 2: The oriented contour 7 is composed by a forward 
and a backward branch between and 00 and a vertical track 
going from to — i/3. According to the orientation the point 
z is later than z' and any point lying on the vertical track is 
later than both z and z' . For any physical time t we have two 
points t± on 7 at the same distance from the origin. In the 
main text we choose the Greek letter r for z on the vertical 
track. 



The basic quantity in NEGF is the Keldysh-Green's 
function 



[Q(z;z') 



Tr 



T ■ 



-i f dzH(z) 



.(*)4(*0} 



Tr[p] 



(4) 

where Tr denotes the trace over all many-body states. In 
Eq. (4) the integral is over the Keldysh contour 7 of Fig. 
2, T is the contour ordering operator and H(z) = H(t) 
for z — t± while H(z) = H for z = t. The fermion 
operators c m (z) — c m , cjj(z') = do not depend on the 
contour variables z and z'\ the reason of the contour ar- 
gument stems from the need of specifying their positions 
in the contour ordering. The Green's function Q obeys 
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an important cyclic relation 37 ' 38 



S(z;0_) = - e -^g(z;-i(3), 



(5) 
(6) 



which is used as boundary condition for the solution of 
the equation of motion 



i—l-H(z)\g(z;z') = 5(z;z')l, (7) 



with 1 the unity matrix. The Kcldysh-Grccn's function 
reduces to the Matsubara-Green's function for z and z' 
on the vertical track of the contour 7. 

In the NEGF formulation of quantum transport the 
non-local Hamiltonian connecting regions L/R to the 
central region C is treated perturbatively to all orders. 
This is possible by rewriting H as £ + V where 



£ = 



£l 
£ c 
£ R 



V = 



V LC 
Vol V C r 
V RC 



(8) 

and by introducing the diagonal Green's function g which 
obeys the equation of motion below 



{i±l-£(z)}g(z;z') = 6(z;z')l. 



(9) 



The diagonal g allows us to convert Eq. (7) into a Dyson 
equation on 7 



}(z;z')=g(z;z')+ [ dz g(z;z)V(z)G(z;z'), 



(10) 



provided g obeys the same cyclic relations (5-6) as Q. We 
wish to stress that Eq. (10) is an exact equation only if 
the contour 7 includes the forward/backward branches 
and the vertical track. As one can see by inspection, the 
removal of the vertical track from 7 is equivalent to set 
V(z = t) = V a = and hence to start from an initial 
Hamiltonian H° = £° , with £° = £(t < 0). 

To calculate the time-dependent total current one 
needs the lesser component of the Keldysh-Green's func- 
tion. Choosing z = t_, z' = t' + in Eq. (10) and using the 
Langreth theorem 39,40 we find 



g<(t-f 



+ 



poo 

g < (t;t')+ / dtg R (t;t)V(t)g<(t;t') 
Jo 

/>oo 

/ dig < (t;t)V(t)g A (i;t') 
Jo 



-H3 



dfg(t;f)V(f)g(f;t'), 



(11) 



where Green's functions with superscript R/A are re- 
tarded/advanced Green's functions. It is possible to show 
that Eq. (11) is equivalent to 6 



g<(t;t') = g K (t;0)g<(0;0)g A (0:t'). 



Equation (12), as opposed to Eq. (11), has a simple phys- 
ical interpretation. The lesser Green's function is com- 
pletely known once we know how to propagate the one- 
electron states in time and how these states are populated 
before the system is driven out of equilibrium. The time- 
evolution is fully described by the retarded or advanced 
Green's functions g R,A j and the initial equilibrium dis- 
tribution is fully described by g < (0; 0) = if(H°), where 
f(uj) = \/{e l3{ - 0J -^ + 1) is the Fermi function. 

Equation (12) is valid for non-interacting electrons and 
for mean-field (Hartree or Hartree-Fock) electrons; it is 
also valid for Kohn-Sham electrons since in TDDFT the 
electron-electron interaction is described by an exchange- 
correlation potential which is local in time. The present 
work deals with non-interacting electrons (Section III and 
IV) and Kohn-Sham electrons (Section V), and hence Eq. 
(12) can always be used. Due to the non- locality in time 
of the correlation self-energy, Eq. (12) is no longer valid 
beyond the Hartree-Fock approximation of many-body 
perturbation theory. In this latter case the exact formula 
for the lesser Green's function (which fully accounts for 
transient and correlation effects) can be found in Ref. 6. 



III. CURRENT FORMULA FROM NEGF 

The total current from region a — L,C, R can be cal- 
culated from the time derivative of the total number of 
particles in a 



I a {t) = -e-Tr a [-ig < (t;t)] 



(13) 



where Tr Q denotes the trace over a complete set of one- 
particle states for region a and e is the electron charge. 
From the equation of motion (7) the change per unit time 
of the lesser Green's function is proportional to the com- 
mutator [g < (t; t),H(t)]. Projecting the Green's function 
onto different subregions 
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Qll Qlc Qlr 
Qcl Gcc Gcr 
Grl Grc Grr 



(14) 



it is straightforward to realize that 



(12) 



I a {t)=2eReTr c [G^ a (t;t)V aC (t)], a = L,R. (15) 

Equation (15) will be our starting point for the calcu- 
lation of the long-time limit of I a {t) for different initial 
conditions and in the presence of bound states. Equa- 
tion (15) is valid for both non-interacting and Kohn- 
Sham electrons. In the latter case G K is the Kohn-Sham 
Green's function and I a {t) is the total current of the 
Kohn-Sham system which, according to the Rungc-Gross 
theorem, 13 is the same as the total current of the real in- 
teracting system. In the remaining of this Section and 
in Section IV we specialize to non-interacting electrons. 
Interacting systems will be considered in Section V. 
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A. The contacting approach 

The contacting approach has been originally intro- 
duced by Caroli et al. in a series of four papers. 1-4 About 
twenty years after it has been combined with NEGF to 
include the effects of time-dependent perturbations and 
short-range correlations in the scattering region. 5 ' 41 In 
the contacting approach regions L and R are in equi- 
librium at the same temperature and chemical poten- 
tial and are disconnected from C. The Hamiltonian for 
negative times is then H° = £° and the one-particle 
eigenstates of the unperturbed (and hence disconnected) 
system are strictly confined in the left region or in the 
right region or in the central scattering region. The equi- 
librium distribution of region C can not be univocally 
fixed since C is initially an isolated subsystem. This is 
equivalent to say that £c(t < 0) = £° c is in principle 
an arbitrary Hamiltonian. We observe that this am- 
biguity make the transient current difficult to interpret. 
To drive a current through the system one exposes the 
electrons to a longitudinal electric field, and at the same 
time switches on the contacts V between C and L/R. In 
a typical experimental set up regions L and R are bulk 
metallic electrodes, see Fig. 1. The dynamical forma- 
tion of dipole layers at the a-electrode interface screens 
the potential drop along region a and the total poten- 
tial turns out to be uniform in the left and right bulks. 
For this reason we model the effects of a time-dependent 
electric field uniform shift of the energy levels, i.e., 
£ a (t) — £ a + U a (t)l a , a = L, R and l a the unity matrix 
in region a. 

Using Eq. (11) the time-dependent current of Eq. (15) 
can be rewritten as 5 

/>OC 

I a (t) = 2eReTr c / di 
Jo 

x 'G R (t;t)E<(t;t) + Gf < (t;t)S*(t;t)] , (16) 

where we have defined the so called embedding self- 
energy 

S Q (z; *') = V Ca (z)g aa (z; z')V aC (z'), (17) 

and introduced the short-hand notation 

G(z;z') = g cc (z;z') (18) 

for the Green's function projected in region C. In obtain- 
ing Eq. (16) we have taken advantage of the fact that 
V(z) vanishes on the vertical track (isolated subsystems) 
and hence the last term in Eq. (11) does not contribute. 
For the same reason the embedding self-energy vanishes if 
z and/or z' lie on the vertical track and the lesser Green's 
function G < appearing in Eq. (16) can then be rewritten 
as (see Eq. (Al) of Ref. 6) 



with S 



-L.R, 



S Q the total self-energy. In most 



G 



I'OC 

: (t;t') = / didFG R (t;t)-E < (i;F)G A (i';t') 
Jo 

+ G R (t;0)<7< c (0;0)G A (0;O, (19) 



works on quantum transport the second term in Eq. (19) 
is neglected (see for instance Eq. (31) of Ref. 41). This 
is possible in the limit t — > oo and/or t' — > oo provided 
the local density of states (LDOS) Dq in C is a smooth 
function. However, Dq is not smooth if bound states are 
present. Bound states render the contacting approach 
ambiguous because the initial equilibrium of region C 
affects the behavior of the system also for t — > oo. More- 
over, as it has been recently pointed out, 11 the time- 
dependent current and one-particle density matrix do not 
reduce to their equilibrium values when the perturbing 
electric field is set to zero. 

In the next Section we describe the partition-free ap- 
proach. We will obtain a general expression for I a (t) 
which depends explicitly on the initial equilibrium. This 
expression also allows us to switch easily between the 
contacting and the partition-free approaches. 



B. The partition-free approach 

The partition-free approach has been introduced by 
Cini 12 about a decade after the works of Caroli et al. 
Here the system is contacted and in equilibrium at a 
unique temperature and chemical potential before an ex- 
ternal electric field is applied. The initial density matrix 
of the central region is then uniquely defined and tran- 
sient phenomena have a direct physical interpretation. 
Substituting Eq. (12) into Eq. (15) and performing the 
multiplication between Green's functions we obtain 

I a (t) =2eReTi c [Q a (t)}, a = L,R (20) 

with 

Q a (t) = E Sc 7 (*;0)C^(0;0)^, a (0;t)V a c(i) 

77 ; — L.R 

+ E ^c-y(t;0)g< c (0;0)g^ a (0;t)V aC (t) 

1=L,R 

+ E ag c (i;0)a< 7 ,(0;0)a A a (0;t)F aC (i) 

j'=L,R 

+ g^ c (t;0)QccMQc a (0;t)v aC (t). (21) 

Equations (20,21) are completely general. In the 
partition-free approach P < (0;0) = if(H°) = if(£° + 
V ) while in the contacting approach P < (0;0) = 
if(H°) = if(£°) with £% arbitrary. In Appendix A 
we prove that in the contacting approach Eqs. (20,21) 
are equivalent to the current formula in Eq. (16), as it 
should. 

In Ref. 6 we have shown that by exposing the electrons 
to a constant (in time) electric field the total current I a (t) 
in Eqs. (20,21) tends to a steady value given by (for, 
e.g., a = L) 

7 (S) = e | g [/(w _^°)_ /(w _^ )]T(a;); (22 ) 
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with T(uj) = Tv c [G r (lo)T l (uj)G a (lu)T r (lu)}, T a (uj) = 
-2Im[£*(w)], and J7~ = lim^^ t/ a (t). We have also 
shown that the steady value is independent of the his- 
tory of the applied bias (memory-loss theorem) and of 
the initial equilibrium distribution of region C (theo- 
rem of equivalence). These results were obtained for a 
smooth LDOS Dc- The presence of bound states requires 
a generalization of the theory. In the next Section we in- 
clude bound states in time-dependent quantum transport 
and investigate how they affect the behavior of I a (t) for 
t — > oo. 



IV. INCLUDING BOUND STATES IN 
QUANTUM TRANSPORT 

For simplicity we consider the case of a sudden switch- 
ing on of the bias U a (t) — 9(t)U^ in region a = L,R (ar- 
bitrary time-dependent biases are considered in the next 
Section). We also specialize the discussion to a hopping 
Hamiltonian V(z) = V°° that is constant on the for- 
ward/backward branch of the contour 7. However, we 
leave open the possibility of having a different constant 
value of V(z) = V a for z on the vertical track (this allows 
us to deal with different approaches at the same time). 

In quantum transport experiments the central device 
is typically connected to macroscopic metallic electrodes. 
The one-particle cigenstates \<pka) of = £ a {t — > 
00) = £° a + U™l a , a = L,R, form a continuum with 
energy e^ a = e a ka + e W a . We define a bound state 
of the whole biased system L + C + R as an eigenstate 
\ip b ) of H°° = H(t -> 00) = £°° + V°° with energy 
£b Wl U Wr. [For one-particle eigenstates of the iso- 
lated regions £ a , a — L,R, we use the Greek letter <j) 
while for one-particle eigenstates of the connected and 
biased system if 00 we use the Greek letter tp.] Below 
we calculate the long-time limit of Q a in the presence of 
bound states. 



A. Asymptotic kernel Q a 

According to Eq. (21), the long-time limit of the ker- 
nel Q a (t) is known once we know the long-time limit of 
eScfrO), Sc 7 (*;0), and 0£ a (O;i)V~ C) G A , a (0;t)V™ c . 
For any t > the retarded Green's function can be writ- 
ten as 

g K (t;0) = -iexp[-iH°°t] = J ^g R (o;) e --*, (23) 

and Q A (0;t) = [0 R (t;O)]t. Exploiting the smoothness of 
the self-energy, we can use the Riemann-Lebesgue theo- 
rem to derive the following asymptotic behaviors 

lim <? R c (i;0) = -iJ2^bc){i> b c\e- l£bt , (24) 

b 



lim g^(t;0) = -i^ b c)(Ac\V^— 

- J ^G R ( £ r 7 )^ 7 |0 fc7 )(0 fc7 |e-^*,(25) 
k 

and 

lim GLMV™ c =iJ2e iEbt \i> b c)(i> b c\V A (s b ), (26) 
t — >oo * — * 

h 

lim a^ a (O;t)V~ c = i<y ra y;e fe S o -*|^ a >(0 fca |V~ c 

k 

+^ e ^*|0 fc7 ,)(0 fc7 ,|y- c G A ( £ r 7 o^( £ r y ) 

k 

+i S g i _ro V^c\M(i> b c\V A (s b ), (27) 
b £ b±y> c 7 ' 

where \ipbc) is the projection of eigenstate |^>(,) onto re- 
gion C* and G R ' A (w) = (w). 

The above asymptotic results are given in terms of 
sum over bound states (discrete part) and sum over 
the continuum of states of £ a , a = L,R (continuum 
part). Inserting them into Eq. (21) we obtain a 
continuum-continuum contribution, , and a discrete- 
discrete contribution, Q^, as well as cross terms with 
discrete-continuum contributions. Taking advantage of 
the smoothness of G R,A (w) for uj e Wl U Wr and ex- 
ploiting the Riemann-Lebesgue theorem it is possible to 
show that all cross terms vanish for t — > 00. 

Let us focus on the continuum-continuum part of Q a 
(it is straightforward to realize that only the first term 
in Eq. (21) can contribute to this part). Expanding 
g < (0: 0) in Matsubara modes 

g<(0;0) = J-Y / e^g M (cu n ), r,^0 (28) 

— lb L ' 

n 

with g M (uJ n ) = l/(w n l - H°) and u n = n + (2n + 
l)w/(—i(3), we can rewrite £? 77 /(0;0) as 

e< 7 ,(o ; o) = + -L^ e ^ g M K) 

^ n 

xF 7C ^ c K)V° C7 ,^ 7 ,K)' ( 29 ) 

with g 77 (o)„) = l/(w„l 7 — £ J ). We observe that in the 
contacting approach V° = (regions L,C,R isolated for 
negative times) and hence the second term in Eq. (29) 
vanishes. In this case, the continuum-continuum part is 

Qi s) - * J ^/(w - £C)G R Mr Q M 

+* / £ E - c/ 7 oo )G R (^)r 7 ^)G A (c)s A (c),(30) 

which does not depend on time. In Appendix B we prove 
that does not depend on V° provided the hopping 
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matrix elements between states in region C and states 
\4>ka)i — L, R, are smooth functions of k. Thus, the 
memory of different initial conditions is washed out by 
the continuum of states and the final result is a steady 
contribution. Below we show that this is not the case for 
the discrete-discrete part . 

All four terms in Eq. (21) contribute to Q^ D) . After 
some algebra one finds the following dynamical kernel 



Qi D) W = 



6,6' 



/6,6#6C>(VVc|£«(£6')e- 



»(£(,- 



,)t 



with 



fl 



6.6' 



(^|/(Jf )|W). 



(32) 



The coefficient f btb > is the matrix element of the equi- 
librium Fermi function f(H ) between bound states b,b' 
of Hit -> oo ) = H°°. Equations (31,32) have been ob- 
tained using the bound-state contributions of the asymp- 
totic behavior in Eqs. (24-27) and the relation 



IVV) = — 7 
S b l 



(33) 



with IV'67) the projection onto region 7 = L, R of the full 
bound state \ip b ). 

In conclusion we have 



(K) 



lim Q a (t) = Q ( a 

t—>oo 



oL D) «). 



(34) 



Equation (34) is completely general and is valid provided 
the Hamiltonian of region a = L, R is uniformly shifted 
by a constant potential = £° a + U^l a ), the con- 
tacting Hamiltonian V{t) = V°° and the Hamiltonian 
for the central region £c(t) — £c are constant in time. 
We observe that in deriving Eq. (34) no statements about 
V a and £° c have been made; in principle, they might be 
different from their asymptotic values V°° and . 



B. Asymptotic current and one-particle density 
matrix 

The total current in the long-time limit is obtained by 
substituting Eq. (34) into Eq. (20) 



lim I a (t)=I^+lV\t), 



(35) 



where 1^ is the steady part coming from and is 

given in Eq. (22). Taking into account that the imagi- 
nary part of vanishes for u = e b , the dynamical 
part I^\t) can be written as 

IM (t) =2eJ2 fb,b>A$ sin [(e b - e v )t] , (36) 

6,6' 



with 



A 



(a) 
6,6' 



Tr r 



|V>6C}(VVC|S£( £ 60 



(37) 



Proceeding along similar lines one can also prove that the 
one-particle density matrix p c (t) — —iG < (t;t) in region 
C is the sum of steady and dynamical contributions 



Hmp c (t) = p ( § ) +p ( g\t), 



(38) 



with 



Pc^ / ^^f^-U^G^T^G^), (39) 



(31) and 



Pc (*) = ^hA^bc){^b'c\e 



i{e b —e b i)t 



(40) 



6,6' 



Equations (32,36,40) and the following discussion (in- 
cluding Section V) are the main results of this work. 

In the contacting approach (V° = 0) H° = £° is a 
block-diagonal matrix and so is f(£°). The coefficients 
fb.b' can then be rewritten as the sum of three terms 



h 



6,6' 



(V>6c|/(^)|VVc) + ^ (^671/(^7)1^7) 
1=L,R 



(41) 



The first term in f b p depends on the initial equilibrium 
of the isolated central region and can not be univocally 
fixed (see discussion in Section III A). To make contact 
with Ref. 11 we insert a complete set of eigenstates 0fc 7 ) 
for region 7 = L, R in the second term of Eq. (41). We 
find 

<V>6 7 |/(£ 7 )|VV7) = £/(4y) 

k 

(^bc\V^\cj> k7 )(cj> k7 \V^ b ,c) 



(e b -e™){e v -s™) 



-I 



k-y) 



d \f^-U-) { f^ T t^ b ' C l (42) 



2ir- 



(e b -uj){e bl -w) 



which agrees with the result of Dhar and Sen (see Eq. 
(55) of Ref. 11). We also observe that f b y of Eq. (41) 
gives rise to oscillatory terms in the total current and one- 
particle density matrix which do not vanish for — 
0, 7 = L,R. Thus, in the unbiased system I a (t) and 
p c (t) oscillate forever. To solve this problem Dhar and 
Sen have proposed a somewhat ad hoc procedure which 
consists in introducing two extra reservoirs with a wide 
band. The bandwidth should be large enough to allow 
for the hybridization of the originally localized bound 
states. Eventually, the coupling to the extra reservoirs is 
removed and the standard equilibrium result is recovered. 
However, out of equilibrium this procedure suffers from a 
serious problem: The non-equilibrium quantities depend 
on how these extra couplings approach zero. 

The partition-free approach is free from all the limita- 
tions described above. The initial equilibrium is unam- 
biguous and £° c is simply given by the projection of the 
physical Hamiltonian onto region C . For the unbiased 
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system the equilibrium Hamiltonian H and the long- 
time limit Hamiltonian H°° are the same. In this case 
Eq. (32) implies 

fb,v = SbM f(e b ), (43) 

and the dynamical contribution to the total current van- 
ishes while the dynamical contribution to the one-particle 
density matrix reduces to the equilibrium contribution of 
bound states 

P ( c ] =Y.f^)\^c){^ b cl (44) 

b 

as it should. Out of equilibrium, both I a (t — > oo) and 
p c (t — > oo) oscillate around some steady value provided 
H°° has more than one bound state solution. The ampli- 
tude of the oscillations are calculable expressions and are 
completely fixed by the original temperature and chemi- 
cal potential. Non-equilibrium results are well defined. 

In Appendix C we also prove that Eqs. (36,40) con- 
serve the total number of particles 

^N c (t) = - e [I L (t) + I R (t)}, (45) 

where 

N c (t) = Tr c {p c (t)} (46) 
is the number of particles in region C. 

V. IMPLICATIONS FOR TDDFT 

One of the main advantage of the partition-free ap- 
proach over the contacting approach is that the former 
can be combined with Time-Dependent Density Func- 
tional Theory (TDDFT). 6-8 In this theory the time- 
dependent density of an interacting system can be calcu- 
lated from a fictitious system of non-interacting electrons 
moving under the influence of the Kohn-Sham (KS) po- 
tential U. The KS potential is the sum of the external 
applied potential u ex t, the Hartree potential wh and the 
exchange correlation potential i> xc . According to the dis- 
cussion of Section III A, the metallic screening keeps the 
bulk electrodes in local equilibrium and we can approx- 
imate U = (w xt + + ^xc) with a spatially constant 
time-dependent shift deep inside region L/R. The one- 
particle Hamiltonian of TDDFT is then given by Eq. (2) 
with £ a (t) =£° a + U a (t)l a , a = L,R, and V{t) = V . 

Let us expose the electrons to a constant (in time) 
electric field and let us assume that the system reaches 
a steady state in the long-time limit. Then, the effec- 
tive potential of the bulk electrodes and the Hamiltonian 
of the central region £c are stationary in the distant 
future. We have shown that the steady-state assump- 
tion is consistent with the TDDFT equation for the to- 
tal current provided the density of states in region C is 
a smooth function. 7 All history and initial-state effects 



are contained in v xc (r,t — ► oo), meaning that two dif- 
ferent time-dependent densities n(r,t) and n'(r,t) may 
give the same total current. In the Adiabatic Local Den- 
sity Approximation (ALDA) there is no memory and the 
exchange-correlation potential depends on the instanta- 
neous density. Hence, v ALDA (r,i — > oo) is completely 
known once we know n(r, t — > oo). The latter can be cal- 
culated self-consistently from Eq. (39). The ALDA pro- 
vides a practical scheme to compute the current. How- 
ever, such a scheme, originally proposed by Lang, 43 is ob- 
viously limited to the ALDA. Moreover, owing to the non 
linearity of the problem there might be multiple steady- 
state solutions, and the absence of a minimum principle 
in out-of-cquilibrium systems makes impossible to say to- 
wards which steady-state the electrons actually evolve. 

Below we show that the steady-state assumption is not 
consistent in the presence of bound states. This result 
opens up the possibility of having self-consistent oscilla- 
tory solutions even for constant biases and may change 
substantially the standard steady-state picture. Indeed, 
oscillations of the effective potential in region C give rise 
to new conductive channels, very much in the spirit of 
what happens in quantum pumps. 

The proof of the above statement proceeds by reductio 
ad absurdum. Let H°° be the steady-state Hamiltonian 
of the fictitious non-interacting system. In Appendix D 
we prove that the total current I a (t — > oo) and one- 
particle density matrix p c (t — ► oo) oscillate like in Eqs. 
(35,38) but with new coefficients fop. An oscillating 
density /current is not consistent with the steady-state 
assumption; the effective potential of TDDFT is a func- 
tional of the density and hence H(t) can not be constant 
in time. We conclude that a steady current is not com- 
patible with the existence of bound states in the steady 
Hamiltonian H°° . 

The new coefficients 

fop = (^|/(ff°M,> (47) 

depend on the history of the bias. Indeed, the state \tp' b ) 
is related to the bound state \ipb) 01 H°° by a unitary 
transformation 



W bL ) 




~e* A ~l L " 




' \4>b L ) ' 


W bC ) 




M c 




H>bc) 


Wm) 




e iA *l R 




_ \lpbR.) _ 



(48) 
with 

= lim f At(U a (t) - C/ Q °°) , a = L, R, (49) 

and Mc some unitary "memory matrix" for region C. 
For constant KS potentials = and Mc = lc an d 
the coefficients fop reduce to those in Eq. (32). 

The above results can be extended to systems of 
interacting electrons and bosons, like, e.g., phonons. 
We consider the system electrons+phonons initially in 
equilibrium and assume that the electron density is v- 
representable. Then, there must exist a local poten- 
tial that reproduces the interacting electron density in 
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a non- interacting system of only electrons. Such a po- 
tential defines uniquely H°. Next, we drive the system 
electrons+phonons out of equilibrium by switching on a 
longitudinal electric field and we ask the question if the 
interacting time-dependent electron density can be re- 
produced in the non-interacting system of only electrons 
moving under the influence of a time-dependent local po- 
tential U. According to the van Leeuwen theorem 16 the 
answer is affirmative provided the electron-phonon in- 
teraction preserves the continuity equation, which is a 
very weak condition. Moreover, the local potential U is 
unique. Again, we can conclude that if the Hamiltonian 
of the fictitious system globally converges to an asymp- 
totic Hamiltonian H°° for t — ► oo, the system can not 
have more than one bound state. 

We also observe that for truly non-interacting electrons 
the effective potential is equal to the external potential 
and does not depend on the density. In this case the 
solution of Appendix D is an exact solution and Eqs. 
(36,40), with fop given by Eq. (47), can be tested 
against history-dependent effects in total currents and 
densities. All history (of the applied bias) and initial- 
state dependence are contained in the coefficients fop. 

Finally, we wish to discuss a rather interesting exam- 
ple. Let us imagine to switch on and then off the bias in 
a system of (i) truly non-interacting electrons and (ii) KS 
electrons. In system (i) there is no self-consistent evolu- 
tion since the Hamiltonian is independent of the density. 
The asymptotic Hamiltonian H°° is (trivially) equal to 
the initial Hamiltonian H° and hence the bound states 
are also eigenstates of H°. Suppose that the bound- 
state energies lie in the continuum when the bias is on. 
Then, switching the bias off would result in a depopu- 
lation of the bound states and the asymptotic density 
matrix p c (t — > oo) would, in general, differ from its 
equilibrium value. This expected result is confirmed by 
the exact solution. The system "remembers" that a bias 
has been switched on through the memory matrix M c 
and the phases A£°, a = L, R, appearing in Eq. (48). 
Equation (48) defines new states ip' b ^ ipb and according 
to Eq. (47) the coefficients fop are no longer given by 
5b,b'f( e b)i as it would be in equilibrium. The situation 
is totally different in system (ii). The time evolution of 
interacting electrons in initial noncquilibrium states has 
been recently investigated within Time-Dependent Cur- 
rent Density Functional Theory 44 - 45 (TDCDFT). It has 
been shown that the inclusion of memory effects in the 
exchange-correlation vector potential leads to a dissipa- 
tive KS dynamics, 47-49 and D'Agosta and Vignale have 
been able to prove 48 that the electron density evolves (un- 
less forbidden by symmetry) towards the ground-state 
density. Taking into account that the density of the 
TDDFT KS system is the same as the density of the 
TDCDFT KS system we conclude that H°° = H° and, 
more important, that the memory matrix Mc = lc an d 
the phases &. a =L,R = (mod 2n), independently of the 
history of the switching on/off process. Indeed, for these 
values of Mc and &. a =L,R the density matrix is constant 



in time (and equal to the ground-state density matrix), 
and hence compatible with the (stationary) ground-state 
KS Hamiltonian. 



VI. CONCLUSIONS 

We have generalized the theory developed in Rcf. 
6 to include bound states in quantum transport. In 
the partition-free approach the electrode-device-electrode 
system is contacted and in equilibrium at a unique tem- 
perature and chemical potential (thermodynamic consis- 
tency). The electrons are exposed to a time-dependent 
electric field for positive times. The external potential is 
local in space and total currents and density can be calcu- 
lated from a fictitious system of non-interacting electrons, 
according to the Runge-Gross theorem. 

We have shown that for truly non-interacting elec- 
trons the biased system does not evolve toward a steady 
regime. The total current and density oscillate due to 
the presence of bound states. The amplitude of the oscil- 
lations depends on the initial temperature and chemical 
potential and on the history of the applied bias. Bound 
state oscillations might provide a probe to unveil the past 
of the system. 

In contrast to the contacting approach, in the 
partition-free approach the initial equilibrium distribu- 
tion of region C is well defined and all quantities reduce 
to their equilibrium value by setting the external poten- 
tial to zero. There is no need of "extra reservoirs" to 
equilibrate spurious bound-state oscillations. 

Our findings might have interesting implications for 
the fictitious system of TDDFT. In Rcf. 7 we have 
shown that the time-dependent current tends to a steady 
value provided the Hamiltonian globally converges to a 
steady Hamiltonian and the density of states in region 
C is smooth. This result is no longer valid in the pres- 
ence of bound states. In Section V we have shown that 
bound electrons in steady state regimes lead to a con- 
tradiction: current and density would oscillate and, as 
a consequence, the effective potential of TDDFT would 
oscillate too. Steady quantities are not compatible with 
the existence of bound states. 

According to the above discussion, the effective po- 
tential of TDDFT might oscillate upon application of a 
constant bias. Wc expect that these oscillations have ex- 
ponentially small amplitude deep inside the electrodes 
and are detectable only close to the molecular device. 
Indeed, for truly non-interacting electrons the ampli- 
tude of the density oscillations is proportional to the 
bound-state wavefunctions. Therefore, the KS potential 
is time-dependent in the device region and tends expo- 
nentially to a constant (in time) deep inside the elec- 
trodes. In this case one can use a recently proposed 
practical scheme 8 to investigate bound-state dynamical 
effects within TDDFT. The scheme is based on the real- 
time propagation of the occupied KS orbitals and we are 
currently working at the implementation of the algorithm 
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(which has been tested in one-dimensional model sys- 
tems with excellent results 8,9 ). As for the case of gate 
voltages in electron pumping, oscillations of the effective 
potential open new conductive channels. Such an effect 
is completely left out in static DFT calculations and its 
possible relevance in molecular transport has yet to be 
discovered. 
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APPENDIX A: EQUIVALENCE BETWEEN 
CURRENT FORMULAS IN THE CONTACTING 
APPROACH 

In the contacting approach regions L,C,R are isolated 
for negative times and hence Q^q(0; 0) = £/py(0; 0) = 0, 
<?< 7 ,(0;0) = Vffrr^O) and S< c (0;0) = g£ c (0;0). 
The kernel Q a (t) simplifies to 

Q a (t) = J2 g c^0)g^(0;0)g^ a (0;t)V aC (t) 

7 

+ 0cc(t;O)0£ c (O;O)0£ a (O;t)V aC (t). (Al) 

Extracting the retarded/advanced component of the 
Keldysh-Green's function from the Dyson equation (10) 
we can express in terms of Q R q — G R , and Q A 

Qca m terms of Qcc — Equation (Al) can then be 
rewritten as 

Q a (t) = jG R (t;i)?;<(i;F)G A (F;t")VZ(i";t) 

; R (t;0) fl £ c (0;0) J G A (o ; mt(t;t) 

J G R (t;t)-E<(t;t), (A2) 



+ G 

+ 



where the integrals (between and oo) are over barred 
time variables. Taking into account Eq. (19) it is 
straightforward to realize that Eq. (16) is actually equiv- 
alent to Eqs. (20,21). 



APPENDIX B: INDEPENDENCE OF V° IN Q^ s) 

The continuum-continuum contribution to Q a (t) can 
be written as + SQ a . The steady value orig- 
inates from the first term on the right hand side of Eq. 



(29) and is independent of V . The extra term SQ a ac- 
counts for the possible effects of a non- vanishing coupling 
V a in the remote past. Here, we prove that SQ a vanishes 
in the long time limit. From Eq. (29) we find 



5Qa = 4gE e ^ n ^aK)' 



(Bl) 



with 



7 J 



dm G r (uj)S^(uj) 
2^uj n - oj + U™ 



e-^*G M K) 



„5 aY +G A (to>)?: A (u')^ t 



The matrix S 7 (cj) is defined according to 



.(B2) 



S» =J22n8{u,-e%)V% \^)(^\V° jC , (B3) 

k 



and is a smooth function of u> provided the matrix ele- 
ments of V° between states in region C and states |0fc 7 ), 
7 = L,R are smooth functions of k. Also, S 7 (w) = for 
lu £ W L UW R and hence the products S' 7 G R and Sl /I G A 
are smooth functions. At finite temperature the Mat- 
subara frequencies have a finite imaginary part and the 
denominators in Eq. (B2) are well behaved. Exploiting 
the Riemann-Lebesgue theorem we conclude that both 
integrals in Eq. (B2) vanish for t — > oo, meaning that 

SQ a = o. 



APPENDIX C: CONSERVATION OF THE 
PARTICLE NUMBER 

The time derivative of the number of particles in region 
C can be easily calculated from Eq. (40) 



dt 



N c (t) = - E( £fc - £ b>)fbM>S b M> sin [(e b - e b >)t] 
b,b> 



(CI) 

with S b ,b' the overlap between bound states b, b' in C 

S h ,v =Tr c [|V&c)<Vvc|]. (C2) 

Let us now consider the expression for the total current. 

fs) (S) 

The sum I L + I R of steady contributions vanishes. 
Taking into account that fop is symmetric under the 
exchange of b and &', and that the real part of the ad- 
vanced self-energy is a Hcrmitian matrix we can safely 
replace A["2 with 



b.b> 
1 



t(e b ))} (C3) 



in the expression for the dynamical contribution to I a (t) . 
The difference between self-energies at different bound- 
state energies is 

(£6 - £6') 



K(sb')-K(e b ) = v 



Ca [e v l a -£Z][e b l a -eZ] aC ' 

(C4) 
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and hence Eq. (C3) becomes 

A (o) = (jb_^d TTa j |^ 6a)( ^, a |] ; (C5 ) 

where we have used Eq. (33). Exploiting the orthonor- 
mality relation 

^[IV^X^'al] =h,v, (C6) 

a=L,C,R 

we find 

E Ag, = -^^S M „ (C7) 

It is now straightforward to realize that the sum 1^ + 

1^ of dynamical contributions is equal to the change per 
unit time of the number of particle in region C [which is 
given in Eq. (CI)]. 



For large t we can replace £c(t) with its asymptotic value 
£q. Moreover, taking into account that the self-energy 
vanishes when the separation between its time arguments 

R — R 

goes to infinity we can replace X^(i;f) with H a (t;t), in 
accordance with Eq. (D3). Then, Eq. (D4) reduces to 
the equation obeyed by G cc {t — * oo; 0), meaning that 

lim g* c (t; 0) = lim S* c (t; 0)M^, (D5) 

where M^ c is some unitary matrix that accounts for the 
history of the applied bias (memory effects). From Eq. 
(D5) and the equation of motion for Gc a (®>t) one can 
also prove that 

lim 9L(Q; t)V° aC = lim McGtaiO; t)V° aC . (D6) 



APPENDIX D: LONG-TIME LIMIT OF TDDFT 

We have already shown in Section IV that for a sud- 
den switching on, H(t > 0) = H°°, Vi > 0, the total 
current and one-particle density matrix oscillate in the 
long-time limit. Let us denote with a bar (g, Q, and S) 
Green's functions and self-energies corresponding to this 

— R — R 

case. The asymptotic expressions of Q cc (t; 0), G Cl (t; 0), 

and Qc a (0;t)V° a c, G* a (0;t)V° a c are given in Eqs. (24- 
27). Below, we will find a relation between the asymp- 
totic behavior of Q and Q, with Q the Green's function 
of the TDDFT Hamiltonian. 

For arbitrary time-dependent potentials in region a = 
L, R, the retarded self-energy can be expressed in terms 
of S R as 

E£(i; t') = e-^C'S^t; t')e lA ~ {t '\ (Dl) 

with 

A a (t)= / d* (U a (t) - U?) . (D2) 
Jo 

The quantity A a (t) depends on the history of the applied 
bias. We assume U a (t) to approach rapidly enough 
that A a (t) converges to some value for t — > oo. 
Then, 

lim £*(*;*') = S*(*;0- (D3) 

t,t'— >oo 

The above identity allows us to fix the asymptotic be- 
havior of the Green's function in region C. From the 
equation of motion we have 

|iAi c -£ c ( t) Jgg c (t;0) 

- JdiY*;(t;t)g* c (i;0) = 8(t). (D4) 



Also the asymptotic behavior of can be calculated 
from the equation of motion. We have 

{i^lc - £c(*)} g^(t;0) - V° C7 g*(t;0) 

-y"dtS£(t;t)0g 7 (i;O)=O. (D7) 

Taking the limit t —* oo and exploiting the relation 
^(i -» oo; 0) = e^ A ~ g* 7 (i -> oo; 0) we obtain 

lim 0« (t; 0) - e-* A ~ lim 0* (t; 0). (D8) 

t^oo ' t— >oo ' 

In a similar way one can prove that 

lim G$ a {0;t)V° aC = e^> lim 0^ a (0;t)V° aC . (D9) 



Substituting these results [Eqs. (D5-D6) and Eqs. 
(D8-D9)] in Eqs. (20-21) one can calculate the asymp- 
totic behavior of the time-dependent total current. Pro- 
ceeding along the same line which leads to Eq. (35) one 
can show that I a (t) is again given by the sum of the 

steady-state value 1^ of Eq. (22) and the dynamical 
contribution I^\t) of Eq. (36). However, the coef- 
ficients fb,b' are in general different from those in Eq. 
(32). After some algebra one readily finds that the new 
coefficients fb,b' can be expressed as in Eqs. (47-48). 

The time-dependent one-particle density matrix also 
has the same analytic form of Eqs. (38-40) but with the 
same new coefficients fb,v of Eqs. (47-48). 
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